# This script produces Figure 3 in the paper
# Comparison of Party Rank order shift -- Mean vs Median
# Directories set and libraries loaded in script 1.
# Last updated 6th September 2018

set.seed(1234)

rm(list=ls())

## Source in functions script
source('Functions_Mode_Agree_Boot.R')

## Load in CHES Data in R --- these are the raw CHES data
load("CHESData_clean.Rdata")

## Create a function that counts rank flips for mean, mode and median
rankflip <- function(data,question,country,partyname){
	
	expnames <- names(by(data[, question],paste(data[,country], data[,partyname], sep="_"),mean,na.rm=T))
	expmean <- as.vector(by(data[, question],paste(data[,country], data[,partyname], sep="_"),mean,na.rm=T))
	expmode <- as.vector(by(data[, question],paste(data[,country], data[,partyname], sep="_"),Mode,na.rm=T))
	expmedian <- as.vector(by(data[, question],paste(data[,country], data[,partyname], sep="_"),median,na.rm=T))
	expmmm <- data.frame(expmean, expmode, expmedian)

	partycountry <- as.data.frame(matrix(unlist(strsplit(expnames,split="_")),ncol=2,byrow=T))
	expmmm$country <- as.vector(as.numeric(partycountry[,1]))
	expmmm$party <-as.vector(partycountry[,2])

	rankflipmat <- matrix(NA,nrow=length(unique(expmmm$country)),ncol=4)
	colnames(rankflipmat) <- c("country","meanmode","meanmedian","medianmode")

	for (i in unique(expmmm$country)){
		
		dta <- expmmm[expmmm$country==i,]
		dta <- dta[order(dta$expmean),]
		samerank.meanmode   <- sum(order(dta$expmean)!=order(dta$expmode))
		samerank.meanmedian <- sum(order(dta$expmean)!=order(dta$expmedian))
		samerank.medianmode <- sum(order(dta$expmedian)!=order(dta$expmode))

		rankflipmat[i,1] <- i
		rankflipmat[i,2] <- samerank.meanmode
		rankflipmat[i,3] <- samerank.meanmedian
		rankflipmat[i,4] <- samerank.medianmode
	}

	b.rankflip <- colSums(rankflipmat[,2:4]>=1)/length(unique(expmmm$country))

	return(b.rankflip)

}

sims <- 500

sim.rankflip.q1.99 <-matrix(NA,nrow=sims,ncol=3)
sim.rankflip.q8a.99 <-matrix(NA,nrow=sims,ncol=3)
sim.rankflip.q8c.99 <-matrix(NA,nrow=sims,ncol=3)

for (i in 1:sims){
	sim.rankflip.q1.99[i,] <- rankflip(dta1999,"q1","country","partyname")
	sim.rankflip.q8a.99[i,] <- rankflip(dta1999,"q8a","country","partyname")
	sim.rankflip.q8c.99[i,] <- rankflip(dta1999,"q8c","country","partyname")
}

min.pos.99 <- min(sim.rankflip.q1.99[,1])
min.lrg.99 <- min(sim.rankflip.q8a.99[,1])
min.gal.99 <- min(sim.rankflip.q8c.99[,1])

sim.rankflip.q1.02 <-matrix(NA,nrow=sims,ncol=3)
sim.rankflip.q14.02 <-matrix(NA,nrow=sims,ncol=3)
sim.rankflip.q16.02 <-matrix(NA,nrow=sims,ncol=3)

for (i in 1:sims){
	sim.rankflip.q1.02[i,] <- rankflip(dta2002,"q1","country","party_na")
	sim.rankflip.q14.02[i,] <- rankflip(dta2002,"q14","country","party_na")
	sim.rankflip.q16.02[i,] <- rankflip(dta2002,"q16","country","party_na")
}

min.pos.02 <- min(sim.rankflip.q1.02[,1])
min.lrg.02 <- min(sim.rankflip.q14.02[,1])
min.gal.02 <- min(sim.rankflip.q16.02[,1])


sim.rankflip.q1.06 <-matrix(NA,nrow=sims,ncol=3)
sim.rankflip.q10.06 <-matrix(NA,nrow=sims,ncol=3)
sim.rankflip.q12.06 <-matrix(NA,nrow=sims,ncol=3)

for (i in 1:sims){
	sim.rankflip.q1.06[i,] <- rankflip(dta2006,"Q1","country_id","party")
	sim.rankflip.q10.06[i,] <- rankflip(dta2006,"Q10","country_id","party")
	sim.rankflip.q12.06[i,] <- rankflip(dta2006,"Q12","country_id","party")
}

min.pos.06 <- min(sim.rankflip.q1.06[,1])
min.lrg.06 <- min(sim.rankflip.q10.06[,1])
min.gal.06 <- min(sim.rankflip.q12.06[,1])


sim.rankflip.pos.10 <-matrix(NA,nrow=sims,ncol=3)
sim.rankflip.lrg.10 <-matrix(NA,nrow=sims,ncol=3)
sim.rankflip.gal.10 <-matrix(NA,nrow=sims,ncol=3)

for (i in 1:sims){
	sim.rankflip.pos.10[i,] <- rankflip(dta2010,"position","country","party_name")
	sim.rankflip.lrg.10[i,] <- rankflip(dta2010,"lrgen","country","party_name")
	sim.rankflip.gal.10[i,] <- rankflip(dta2010,"galtan","country","party_name")
}

min.pos.10 <- min(sim.rankflip.pos.10[,1])
min.lrg.10 <- min(sim.rankflip.lrg.10[,1])
min.gal.10 <- min(sim.rankflip.gal.10[,1])


pos <- c(min.pos.99,min.pos.02,min.pos.06,min.pos.10)
lrg <- c(min.lrg.99,min.lrg.02,min.lrg.06,min.lrg.10)
gal <- c(min.gal.99,min.gal.02,min.gal.06,min.gal.10)

barplotdta <- t(as.data.frame(rbind(pos,lrg,gal)))
rownames(barplotdta) <- c("Y1999","Y2002","Y2006","Y2010")
colnames(barplotdta) <- c("EU Position","Left-Right","GAL/TAN")

pdf(file = "chesrankorder.pdf", width=12,height=6)
par(mfrow=c(1,2))
barplot(barplotdta, beside=T, xlab="Dimension",legend = rownames(barplotdta), args.legend=list(x="topleft"), main="Percent Countries with Party Rank Order Shift\n Mean vs. Mode")


########################
# Mean Median Comparison

pos.99  <- rankflip(dta1999,"q1","country","partyname")[2]
lrg.99  <- rankflip(dta1999,"q8a","country","partyname")[2]
gal.99  <- rankflip(dta1999,"q8c","country","partyname")[2]

pos.02  <- rankflip(dta2002,"q1","country","party_na")[2]
lrg.02  <- rankflip(dta2002,"q14","country","party_na")[2]
gal.02  <- rankflip(dta2002,"q16","country","party_na")[2]

pos.06  <- rankflip(dta2006,"Q1","country_id","party")[2]
lrg.06  <- rankflip(dta2006,"Q10","country_id","party")[2]
gal.06  <- rankflip(dta2006,"Q12","country_id","party")[2]

pos.10 <- rankflip(dta2010,"position","country","party_name")[2]
lrg.10 <- rankflip(dta2010,"lrgen","country","party_name")[2]
gal.10 <- rankflip(dta2010,"galtan","country","party_name")[2]

pos <- c(pos.99,pos.02,pos.06,pos.10)
lrg <- c(lrg.99,lrg.02,lrg.06,lrg.10)
gal <- c(gal.99,gal.02,gal.06,gal.10)

barplotdta <- t(as.data.frame(rbind(pos,lrg,gal)))
rownames(barplotdta) <- c("Y1999","Y2002","Y2006","Y2010")
colnames(barplotdta) <- c("EU Position","Left-Right","GAL/TAN")

barplot(barplotdta, beside=T, xlab="Dimension",legend = rownames(barplotdta), args.legend=list(x="topleft"), main="Percent Countries with Party Rank Order Shift\n Mean vs. Median")
dev.off()
